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A method and apparatus for accurately and precisely determining the position coordinates of a receiver (15) located 
at a predetermined site using signals broadcast from a plurality of orbiting satellites (13) and using information derived 
from a second receiver (17) located at a reference site having known coordinates. Carrier phase measurements made by the 
two receivers on all of the incoming satellite signals are incorporated into a special square root information filter matrix 
(3 1) that permits a convenient determination of phase noise in the satellite and receiver clocks, for use in subsequently cor- 
recting the phase measurements. The position coordinates of the unknown site are determined in a special three-stage 
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METHOD AND APPARATUS FOR PRECISION 
SURVEYING USING BROADCAST SATELLITE SIGNALS 

Background of the Invention 
This invention relates generally to surveying using signals 
broadcast from a plurality of orbiting satellites, and, more 
particularly, to satellite-based surveying systems that 
determine the position coordinates of an unknown site 
relative to a reference site whose position coordinates are 
known. 



Satellite-based positioning systems such as the Global 
Positioning System (GPS) are now a highly popular means of 
accurately and precisely determining the position of a 
receiver. These systems have numerous practical 
applications and, depending on the time duration over which 
15 measurements are taken, they can determine a receiver's 
position to sub-centimeter accuracy. 

. In GPS, a number of 'satellites orbiting the earth in well 
defined polar orbits continuously broadcast signals 
indicating their precise orbital positions. The broadcast 
signals all have a common- frequency, but are modulated by 
unique, pseudorandom digital codes. Each satellite signal 
is based on a precision internal clock. The receivers 
detect the superimposed modulated carrier signals and 
determine either or both of the code phase and carrier phase 
25 of each detected signal, relative to their own internal 
clocks. These detected phases can be used to determine the 
receivers' position coordinates. 



One typical system for processing these carrier phase 
measurements is described in an article by Bossier et al, 
entitled "Using the Global Positioning System (GPS) for 
Geodetec Positioning", Bulletin Geodesioue . Vol. 54, No. 2, 
1980, pages 553-563. This article describes a processing 
technique known as double differencing. In this technique, 
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the carrier phase measurements are collected at each of two 
receivers for a plurality of GPS satellites. One receiver 
is located at a site whose position coordinates are to be 
determined, and the other is located at a reference site 

5 "whose position coordinates are known* In an initial step, 
the carrier phase measurements for each satellite are 
differenced across the two sites. This eliminates by 
cancellation any clock error in the satellite, since the 
effects of such an error would be identical at each site, 

10 Thereafter, one satellite is chosen as a reference 
satellite, and the single difference measurement obtained 
for it is subtracted from the single difference measurements 
obtained for all of the other satellites. This second 
differencing eliminates by cancellation any clock errors in 

15 the receivers at the two sites, since the effects of such 
clock errors would be identical in all of the single 
differences measurement. 

Several additional steps involving the double difference 
measurements are required. First, rough estimates for the 

20 unknown site's coordinates are provided to the apparatus. 
These estimated coordinates are compared with the specific 
satellite orbits obtained from the data broadcast by the 
satellite or from any other suitable source (e.g., the 
National Geodetic Survey) . The expected value of the double 

25 difference measurements can thereby be determined. 
Equations defining the error in the expected double 
difference values are then formed by differencing them from 
the corresponding actual double difference measurement. The 
sensitivity of these error equations to changes in the 

30 estimated coordinates is also determined, whereby the error 
equations can be solved. The equations are solved usually 
in an iterative, least mean square error procedure, to 
determine the particular site coordinates that will minimize 
the total mean square error in the computed double 
difference values. 
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The carrier phase measurements do not reflect the initial 
number of integral carrier cycles present in each 
satellite/receiver link. The double difference 

measurements, therefore, are biased by an unknown integral 
5 number of cycles. The iterative, least squares procedure 
described above must, therefore, also solve for this bias in 
each double difference measurement. 

When the iterative procedure described above has been 
iterated to the point where the position coordinates and 

10 double difference bias values have converged to fixed 
values, the procedure advances to an additional stage in 
which the computed double differences values are adjusted by 
the nearest whole value to the bias solution and then that 
bias state is dropped from the least squares computation. 

15 Typically, the uncertainty in the bias determination is also 
computed and it is dropped from the least squares 
computation only when the uncertainty is acceptably small. 
Otherwise, the iteration is continued. 



20 



The most accurate 'determination of the unknown site's 
position coordinates is obtained when an account is made for 
the correlation between the double difference measurements 
that arises from differencing with respect to the same 
satellite/receiver link. Accounting for this correlation 
significantly increases both the complexity and time 
25 required to complete the computation. In addition, these 
double differencing techniques are not ordinarily effective 
when the signal for the reference site and/or reference 
satellite is lost. 



30 



There is a need for a technique and related apparatus for 
more effectively utilizing the carrier phase measurements 
obtained in a surveying system utilizing signals broadcast 
from a plurality of orbiting satellites and detected by 
receivers located at two or more sites, one having known 
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position coordinates and the other having coordinates to be 
determined. The present invention fulfills this need. 

Summary of the Invention 
The present invention is embodied in an apparatus and 
related method for determining the position coordinates of a 
first predetermined site using signals broadcast from a 
plurality of orbiting satellites in a significantly more 
efficient and accurate manner than previously performed. A 
first receiver is located at the first predetermined site 
and a second receiver is located at a second, reference site 
having known position coordinates. The two receivers are 
adapted to receive the signals broadcast from the satellites 
and measure the carrier phase of each signal relative to 
their respective internal clocks, thus producing a plurality 
of measured carrier phase values. 

These measured phase values are made at each of a succession 
of time points. To eliminate the effect of an unknown 
integral number of whole carrier cycles in each 
satellite/receiver link, the measured phase values are 
differenced over successive time points. This produces a 
plurality of measured phase difference values, which 
indicate the actual range changes that occur across 
successive time points. The apparatus further calculates 
the expected carrier phase of each satellite signal received 
25 by the respective first and second receivers, and 
differences these expected phase values over successive time 
points, to produce a plurality of expected phase difference 
values. Error values are then generated, equal to the 
difference between the expected phase difference values and 
30 the actual measured phase difference values for each of the 
successive time points. 

Each of these error values has a certain sensitivity to 
incremental changes in the position coordinates of the 
associated site and in the internal clocks of the associated 
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receiver and satellite . The equations for the error values 
are expressed as functions of these latter parameters and 
incorporated into a predetermined matrix, at each successive 
time point. In accordance with the invention, the updating 
of the matrix at each time point includes a preliminary step 
of eliminating the dependence of matrix entries on previous 
values of the phase noise parameters for the receiver and 
satellite clocks. This elimination limits growth of the 
matrix, which otherwise would occur, and thus greatly 
improves the apparatus' efficiency. In the last step of the 
technique, the apparatus solves the matrix to determine the 
particular position coordinates for the first site that 
minimize the mean square value of the plurality of error 
values reflected in the matrix. This solution can be termed 
15 a hyperbolic solution, because it represents the best 
estimate of the single intersection point of a plurality of 
hyperbolic surfaces . 



10 



In one optional aspect of the invention, the matrix is 
formed using a triangular Householder algorithm. The matrix 
entries indicating sensitivity to the receiver and satellite 
clocks are located in a predetermined section of the matrix, 
such that when the matrix entries are updated, dependence on 
previously values of the phase noise parameters of the 
satellite and receiver clocks is eliminated simply by 
zeroing out the predetermined section. The matrix thereby 
remains the same size regardless of the number of times it 
is updated. This substantially simplifies the solving of 
the matrix to determine the first site's position 
coordinates . 



The apparatus and procedure described above can represent 
merely the first stage in the determination of the first 
site's position coordinates. In an optional second stage, 
which is similar in many respects to the first stage, the 
process is modified to take into account the -0.5 
correlation between the error values for successive time 
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points. This correlation arises because one of the terms in 
the respective equations defining the error value for 
successive time points is, by definition, the same. Thus, 
the second stage of the process appropriately whitens the 

5 successive error equations, making them substantially 

uncorrelated with each other, prior to incorporation into 
the matrix. Subsequently solving the matrix yields position 
coordinates that are even more accurate than the first stage' 
solution. This second stage solution can be termed a 

10 Doppler pseudorange solution. 

In an alternative embodiment of the invention, the second 
stage of processing can be used alone, without need for the 
first stage. Unless a fairly accurate estimation of the 
first site's position coordinates is initially made, 
15 however, this second stage will likely need to be iterated. 



In an optional, final stage of the procedure, any ambiguity 
in the number of whole carrier cycles in the measured 
carrier phase signals is resolved. This third stage 
utilizes substantially the same algorithm as the second 

20 stage, described above. The only difference is the first 
measured phase signal is scaled by the square root of 2 
before being incorporated into the matrix. The subsequent 
integrated Doppler measurements incorporated into the matrix 
are made to be uncorrelated with each other, i.e., whitened, 

25 as was done in the second stage, except that a 
cross-correlation coefficient of -.707 is used, rather than 
-0.5. Solving the matrix after the information has been 
fully loaded into it provides a solution equivalent to one 
determined with knowledge of the actual number of whole 

30 carrier cycles in each link. Each satellite/receiver link 
is handled in a sequential fashion, with the correct 
resolution of each improving the probability of correctly 
resolving the next. The final solution can be termed a 
resolved lane solution. 
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Other aspects and advantages of the present invention will 
become apparent from the following description of the 
preferred embodiment, taken in conjunction with the 
accompanying drawings, which illustrate, by way of example, 
5 the principles of the invention. 

Brief Description of the Drawings 
Figure 1 is a schematic diagram (not to scale) of a* 
surveying system having two receivers, one located at an 
unknown site and the other at a reference site having known 
10 coordinates, the receivers detecting signals broadcast from 
four orbiting satellites , to determine the coordinates of 
the unknown site; 

Figure 2 is a simplified flowchart showing the operational 
steps performed by the apparatus of the invention in 
15 accurately determining the position coordinates of the 
unknown site of Figure 1; 

Figure 3 is a schematic diagram of a square root information 
filter matrix used by the apparatus of the invention; and 

Figures 4A and 4B are a flowchart depicting in more detail 
20 the operational steps performed in each of the first, second 
and third stages of th flowchart of Figure 2. 

Description of the Preferred Embodiment 
As shown in the accompanying drawings, this invention is 
embodied in an apparatus for accurately determining the 

25 position coordinates of a predetermined, unknown site 11 
using signals broadcast from a plurality of orbiting 
satellites 13. The apparatus is particularly useful as part 
of the Global Positioning System (GPS) . A first receiver 15 
is located at the unknown site 11, and a second receiver 17 

30 is located at a reference site 19 having known coordinates, 
which can be spaced many kilometers from the unknown site. 
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Eight separate links 21 are thereby formed between the four 
satellites and two receivers. 



The four satellites 13 all broadcast carrier signals having 
the same nominal frequency but each modulated by a unique 
pseudorandom digital code. The respective first and second 
receivers 15 and 17 include antennas 23 and 25 for receiving 
the superimposed, incoming modulated carrier signals and the 
receivers measure the phase angle of each. The phase angle 
measurements are transmitted on lines 27 and 29 from the 
respective receivers 15 and 17 to a data processor 31, for 
position determination. 



The phase measurements are made periodically over an 
extended time (e.g., 30 minutes), and can be compacted in a 
conventional fashion. For example, the measurements can be 
made every 20 milliseconds and compacted into average values 
that are updated once per minute. Phase noise in the 
internal clocks of both the satellites and the receivers 
introduce uncertainties into these successive phase 
measurements . 

In the past, the carrier phase measurements for each of the 
satellite/receiver links 21 were typically processed using a 
double differencing technique to detect the unknown site's 
position coordinates. This technique cancelled out the 
corrupting effects of phase noise in the internal clocks of 
the respective satellites and receivers and thereby led to a 
reasonably accurate position determination. For several 
reasons, however, this prior double differencing technique 
has not proven to be entirely satisfactory. 

The apparatus of the invention avoids the drawbacks 
associated with the prior double differencing technique by 
expressly solving for the magnitude of the phase errors in 
the internal clocks of the satellites 13 and the receivers 
15 and 17, rather than differencing those phase errors out. 
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This yields substantially improved results, in a highly 
efficient manner. 



More particularly, the data processor 31 processes the 
successive sets of measured phase signals in a special 

5 three-stage process, which is depicted schematically in 
Figure 2. Three successive stages 33, 35 and 37 of the 
technique each yield progressively more accurate position 
determinations. The first stage yields what is called a 
hyperbolic solution, the second stage a Doppler pseudorange 

10 solution, and the third stage a resolved lane solution. 
These three stages are depicted in greater detail in Figures 
4A and 4B. The three stages are substantially similar to 
each other, differing in only a few small, but significant 
ways, as will become apparent from the following 

15 description. 

All three of the successive stages 33, 35 and 37 of the 
process expressly solve for the phase errors contributed by 
noise in the internal' clocks of the satellites 13 and the 
receivers 15 and 17.' Each stage determines the difference 
between the phase measurements of each satellite/receiver 
link 21 at successive time points, which corresponds to the 
range change between such time points. Each stage further 
computes the expected range change between the same time 
points, based on the best current estimate of the receiver 
25 site's coordinates and the known satellite orbits. The 
difference between corresponding ones of these phase 
difference measurements and phase difference estimates 
represents error values, whose equivalent equations are 
incorporated into a special square root information filter 
30 (SRIF) matrix. Each equation for each error value is given 
as follows: 



20 



Error (E) = Measured Range - Computed Range 



9E AClock ^ 3e Ann<-v 

3Clock sat «t + 3 clock_.„ lQCk rcvr 



rcvr 
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Where 




Clock 



Clock 



X, Y + Z = Site Coordinates 



= Satellite clock phase noise 
= Receiver clock phase noise 



This equation is formed for all eight satellite/receiver 
links 21. The eight equations include 12 parameters, 
namely, the X, Y and Z coordinates for both receivers 15 and 
17, the phase noise of the four satellite clocks, and the 



reference receiver 17 has known coordinates and is defined 
to have a fixed reference clock, only eight of these 
parameters are unknown. They can be determined by solving 
the eight simultaneous equations. The bias terms in the 
eight equations are cancelled out by subtracting successive 
equations. 

The SRIF matrix, which is depicted in simplified form in 
Figure 3, greatly simplifies the mathematical steps required 
to properly process the information. The lower-left 
triangular section 39 of the matrix is cleared to zero and 
remains in that state for the entire procedure. In each of 
the three stages of the Figure 2 process, the error 
equations for each of the successive time points are 
incorporated sequentially into the upper-right, triangular 
section 41 of the matrix. The particular entries that 
reflect sensitivity of the error values to phase noise in 
the internal clocks of the satellites 13 and the receivers 
15 and 17 are intentionally located in the top rows of the 
matrix. The bottom rows reflect the sensitivity of the 
error values to the site coordinates. By so arranging the 
matrix entries, the dependence of the entries on previous 
values of the clock phase noise can be conveniently 
eliminated simply by zeroing out the pertinent top rows. 
Because of the nature of the SRIF matrix, this has no effect 
on the continued accuracy of the lower rows. In the last 
step in each stage of the procedure, the SRIF matrix is 
solved, to produce the particular corrections to the 



phase noise of the two receiver clocks. 



Because 



the 
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previous estimate of the site coordinates that provide the 
minimum mean square error for the sum of the various error 
values . 



10 



15 



20 



25 



30 



In the first stage 33 of the procedure, the error equations 
are incorporated directly into the SRIF matrix. The 
resulting hyperbolic solution s a very accurate estimation 
of the known site's coordinates. In the second stage 35, 
the equations are whitened (i.e., made to be uncorrelated 
with each other) prior to incorporation into the matrix. 
This yields a slightly more accurate, Doppler pseudorange 
solution. Finally, in the third stage 37, the previous 
uncertainty in the number of whole carrier cycles n each 
satellite/receiver link 21 is expressly resolved, in a 
sequential fashion. The -resulting final solution to the 
SRIF matrix is thereby an extremely accurate estimation of 
the unknown site's coordinates. 

Figures 4A and 4B together are a simplified flowchart 
showing the operational steps performed by the date 
processor 31 in processing the successive carrier phase 
measurements derived for the various satellite/receiver 
links 21. The three stages 33, 35 and 37 of the overall 
process are depicted together in one figure so that the 
extensive similarities of the three stages can be readily 
determined. These similarities lead to significant savings 
in software development and in the sizing of the required 
processor. 

Addressing first the steps performed in implementing the 
first stage 33 of the technique, it will be observed that a 
first step 43 of the procedure clears all of the locations 
in the SRIF matrix. Thereafter, at step 45, the matrix 
locations corresponding to the coordinates of the reference 
site 11 (Figure 1) are constrained to their known values by 
entering appropriate large values, e.g., one million. These 
matrix locations are located along the main diagonal of the 
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SRIF matrix, as indicated by the reference numeral 47 in 
Figure 3. The remaining matrix locations in the same rows 
retain their zero values. 

In a following step 49 of the first stage 33, the program 
5 indexes to accept carrier phase measurements corresponding 
to the first of the plurality of time points. Subsequently 
in the program, when a loopback is made to this step 49, the 
program indexes to the next succeeding time point. 

Next, in a step 51, the processor 31 clears the top rows o 
10 the SRIF matrix that carry information relating to the phase 
noise contributed by the various satellite and receiver 
clocks. In the initial loop through the program, these rows 
are already set to zero (by step 43). Also in step 51, the 
matrix location corresponding to phase noise contributed by 
15 the reference site clock is constrained by entering a large 
value (e.g., one million) into it. This location is the 
diagonal element identified by the reference numeral 53 in 
Figure 3 . 

The program next proceeds to a step 55, where it indexes to 
20 accept the particular carrier phase measurements 
corresponding to the first site 11. Later in the program, 
when looping back to this step, it indexes to the next 
succeeding site. In the preferred embodiment, just two 
sites are described, but it will be appreciated that the 
25 technique of the invention also has utility in processing 
the carrier phase measurements for a greater number of 
sites. 

Next, at a step 57, the program indexes to accept carrier 
phase measurements corresponding to the first of the 
30 satellites 13. Later in the program, when looping back to 
this same step 57, the program indexes to the next 
succeeding satellite. In the preferred embodiment, four 
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13 

satellites are depicted, but it will be appreciated that any 
plurality will suffice. 

At a succeeding step 59, the program computes the range from 
the currently-indexed satellite 13 to the currently-indexed 
receiver site 11 or 19, along with derivatives of the range 
equation with respect to the site's currently-estimated 
position coordinates. Initially, these estimated 

coordinates need not be accurate at all. In fact, they can 
even be initially set at the center of the earth. When the 
estimate is inaccurate to that degree, however, iteration of 
the first stage 33 of the technique will likely be required. 
Preferably, the initial estimate of the receiver's position 
is provided using the codes present on the incoming 
satellite signals. When this is done, iteration of the 
15 first stage is not required. 

The derivatives computed in the steps 59 are used as the 
coefficients in a conventional Taylor series expansion of 
the range equation. First order terms are all that are 
ordinarily required, so that the range equation is a mere 
20 linear function of the estimated site coordinates. 

In a following step 61, it is determined whether or not the 
carrier phase measurements currently being processed 
correspond to the first time point. If they do, the program 
proceeds to a step 63, where it is determined whether or not 

25 the measurements for all of the satellites 13 have been 
processed for the currently-indexed site. If they have not 
been, the program loops back to the step 57, where it 
indexes to the next succeeding satellite. If, on the other 
hand, the measurements for all of the satellites have been 

30 processed, the program advances to a step 65, where it is 
determined whether or not the measurements for all of the 
sites have processed. If they have not, the program loops 
back to the step 55, where it indexes, to the next succeeding 
site. If, on the other hand, the measurements for all of 
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the sites have been processed, the program loops back to the 
step 49, where it indexes to the next succeeding (i.e., the 
second) time point. 

With reference again to step 61, if it is determined that 
5 the phase measurements currently being processed are for the 
second and subsequent time points, the processor 31, at step 
67, computes the differences since the previous time points 
of the phase measurements, computed ranges, and derivatives, 
for the satellite/receiver site currently being processed, 
10 These differences are used to form the error term and 
measurement equation. 

In a following step. 69, a triangular Householder algorithm 
is used to incorporate the measurement equation into the 
SRIF matrix. This is a conventional algorithm that is 

15 described in detail in a book by Gerald J. Bierman, entitled 
"Factorization Method for Discrete Sequential Estimation", 
Academic Press, 1977. The algorithm Bierman describes is 
preferably modified ' to permit operation on a single 
measurement equation at a time, rather than all 

20 simultaneously. 

After the SRIF matrix has been incorporated with the 
measurement equation in step 69, the program advances to 
step 71, where it is determined whether or not the measured 
phase signals for all of the satellites 13 (for the current 

25 site 11 or 19 and current time point) have been processed. 
If not, a return is made to the step 57, where the program 
indexes to the next succeeding satellite. If, on the other 
hand, the measured phase signals for all of the satellites 
have been processed, the program advances to a step 73, 

30 where it is determined whether or not the signals 
corresponding to all of the receiver sites (for the current 
time point) have been processed. If not, a return is made 
to step 55, where the program indexes to the next succeeding 
site. 
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If, on the other hand, it is determined at step 73 that the 
signals corresponding to all of the receiver sites have been 
processed, the SRIF matrix is solved at step 75, to 
determine the phase errors of the various satellite and 
5 receiver clocks. The measured phase signals for the current 
time point are then all corrected to reflect these clock 
phase errors . 

Thereafter, at step 77, the program determines whether or 
not the measured phase signals for all of the successive 
10 time points have been processed. If not, a return is made 
to step 49, where the program inde::es to the next succeeding 
time point. 

If, on the other hand, it -is determined at step 77 that the 
measured phase signals for all of the successive time points 

15 have been processed, the program advances to a step 79, 
where the SRIF matrix is solved to determine the appropriate 
corrections to the site coordinates. These corrections are 
then made to the earlier best estimate of these coordinates, 
and the first stage' of the procedure has been completed. 

20 The resulting corrected position coordinates for each 
unknown site, which, in the preferred embodiment, is merely 
the single site 11, in fact corresponds to the particular 
coordinates that represent the least mean square error in 
the plurality of error signals computed for the various 

25 measurements. These coordinates represent the particular 
location that is the best (i.e., least mean square error) 
estimate of the intersection between the various 
hyperboloids defined by the range difference equations. 

The corrected position coordinates for the unknown site 11 
provided by the first stage 33 of the procedure are 
30 sufficiently accurate for many practical applications. 
Accuracy to within about 5-10 cm can be expected from data 
accumulated over an approximately 30 minute span. One 
factor contributing to the limited accuracy is that the 



SUBSTITUTE SHEET 



WO 87/06410 PCT/US87/00825 



10 



20 



30 



16 

successive time difference or Doppler range equations are 
not entirely uncorrelated with each other. This is because 
each successive difference measurement has, as one of its 
terms, the range equation for the same satellite, receiver 
and time point. The correlation is -0.5. The corrections 
to the position coordinates determined by solving the SRIF 
matrix could be improved if this correlation were 
eliminated. This is what is accomplished in the second 
stage 35 of the procedure. 



As previously mentioned, the second stage 35, which yields 
the Doppler pseudorange solution, is very much similar to 
the first stage 33, which yields the hyperbolic solution. 
With reference to Figures 4A and 4B, it will be observed 
that the second stage differs from the first stage in only 
15 two minor, but significant, respects. First, an extra step 
81 is included immediately after the step 45 of constraining 
the coordinates of the reference site. In particular, this 
extra step clears the correlation coefficient for each 
satellite/site pair.. This is done so that the successive 
difference measurements can be appropriately adjusted to 
make them substantially uncorrelated with each other. 



The second difference between the second stage 35 and the 
first stage 33 is the inclusion of two additional steps 83 
and 85 immediately after the formation at step 67 of the 
25 error terra and measurement equation. In particular, step 83 
incorporates the clock phase noise portion of the current 
measurement into the SRIF matrix. Step 85 then whitens the 
measurement by computing the suitable measurement weighting 
factor from the correlation coefficient and adjusting the 
measurement by combining it with the previous measurement 
multiplied by the correlation coefficient. The whitening 
procedure is described in the Bierman reference identified 
above. The measurement is also scaled by the weighting 
factor. The adjusted and scaled measurement is saved for 
use at the next time point. After step 85, stage two 35 of 
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the procedure continues with step 69 and subsequent steps, 
in the same manner as was described above with respect to 
stage one 33. 

Solving the SRIF matrix at step 79 of stage two 35 provides 
5 corrections to the position coordinates of the unknown site 
11 to make the coordinate estimate even more accurate than 
after stage one 33. Again , this increased accuracy results 
from taking into account the -0.5 correlation between the 
successive signal differences for each satellite/receiver 
10 link 21. 

Although the corrected site coordinates resulting from the 
data processor's implementation of stage two 35 of the 
procedure are sufficiently accurate for many practical 
applications, the accuracy of these coordinates can be 

15 enhanced even further. This enhancement can be achieved by 
taking into account the whole cycle bias in each of the 
successive measured phase signals for all of the 
satellite/receiver links 21. It will be recalled that these 
biases were automatically cancelled out in both stage one 33 

20 and stage two 35 by subtracting the error signals for 
successive time points for each link. Although the biases 
were cancelled out, this particular technique for doing so 
led to minor inaccuracies in the resulting coordinate 
corrections. Stage three 37 of the procedure specifically 

25 solves for the bias in each link. 

With reference again to Figures 4A and 4B, it will be 
observed that stage three 37 of the procedure is very much 
similar to stages one and two, described above. These 
differences are described below. 

30 First, stage three 37 includes an additional step 87 that 
immediately follows the stage two step 81 of clearing the 
correlation coefficient for each satellite/receiver pair. 
In particular, step 87 provides an additional indexing 
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through the plurality of satellites. The first time the 
step is implemented, the program indexes to the first 
satellite. Thereafter, when looping back to this step, the 
program indexes to the next succeeding satellite . Following 
5 step 87, the program continues with steps 49-61, as 
described above with respect to stage one 33, 

In a further difference between the third stage 37 and the 
earlier-described second stage 35, additional steps 89 and 
91 are included following a no answer in step 61, i.e, a 

10 determination that the processor 31 is not currently 
processing signals for the first time point. At step 89, it 
is determined whether or not the satellite number in the 
extra loop (selected in step 8 7) is less than the satellite 
number selected in the mairi loop (in step 57) . If it is, a 

15 return is made to step 49, where an index is made to the 
next time point. On the other hand, if it is determined 
that the extra loop is not less than the first loop, the 
program proceeds to step 91, where the error term is 
computed by keeping 'only the fractional portion of the 

20 difference between 'the measured value signal and the 
computed range. The measurement equation is then formed and 
the correlation coefficient set to -0.707. Thereafter, the 
program proceeds to steps 83 and 85, which is the same as in 
the second stage 35 of the procedure. 

25 The program then continues in the same fashion as in the 
second stage 35, resulting ultimately in the solving of the 
SRIF matrix at step 79. In this third stage 37, however, 
this last step 79 is followed by a step 93, in which it is 
determined whether or not the extra loop through the 

30 plurality of satellites 31 has been completed. If not, the 
program loops back to the step 87, in which the satellite 
number in the extra loop is indexed by one. 

When it is determined at step 9 3 that all of the satellite 
numbers have been indexed through in this extra satellite 
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loop, the entire procedure is completed. The last 
corrections made to the position coordinates of the unknown 
site are, at this time, in their most accurate state. 

The third stage 37 or resolved lane solution is thus divided 
into as many sub-solutions as there are satellites 31. The 
result of each sub-solution is an improved solution that has 
a higher probability of correctly resolving the whole cycle 
ambiguity of the next satellite to be incorporated into the 
solution. The final solution is attained when all of the 
satellites, have been incorporated into the solution and the 
whole cycle ambiguity values resolved. 

Attached as an Appendix is a program listing, written in 
Fortran, of one suitable computer program useful in 
implementing the present invention. 

15 In an alternative embodiment of the invention, the SRIF 

' T 

matrix is replaced by a A A matrix, where A is a standard 
matrix comprising the coefficients of the various error 
equations. This embodiment is less desirable, however, 
because it requires a relatively complicated Gauss 
elimination procedure to eliminate its dependence on prior 
clock phase noise parameters. 

It should be appreciated from the foregoing description that 
the present invention provides a significantly improved 
technique for determining the position coordinates of an 
unknown site using signals broadcast from a plurality of 
orbiting satellites and using the known position coordinates 
of a predetermined reference site. The technique improves 
upon prior techniques by expressly solving for the phase 
noise in the internal clocks of the satellites and 
receivers. The procedure can be broken down into three 
distinct stages, each differing in only a small, but 
significant, way from the preceding stage. Each successive 
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stage yields a slightly more accurate estimation of the 
unknown site's position coordinates. 

Although the invention has been described in detail with 
reference to the presently preferred embodiment , those of 
ordinary skill in the art will appreciate that various 
modifications can be made without departing from the 
invention. Accordingly, the invention is defined only by 
the following claims. 
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I CLAIM: 

1 1. A method for determining the position coordinates of a 

2 first receiver using signals broadcast from a plurality 

3 of orbiting satellites, wherein a second receiver is 

4 located at a reference location having known position 

5 coordinates, and wherein the first and second receivers 

6 each receive the plurality of broadcast signals and 
detect the carrier phase of each signal relative to an 

8 internal clock, the method comprising steps of: 

9 measuring the phase of the cr.rrier signal received by 

10 each of the first and second receivers from each 

11 satellite, to produce a plurality of carrier phase 

12 measurements; 

13 calculating the expected phase of the carrier signal 

14 received by each of the first and second receivers from 

15 each satellite, based on the satellite's orbit and an 

16 estimate of the' receiver's position, to produce a 

17 plurality of carrier phase estimates; 

18 repeating the steps of measuring and calculating at 

19 regular time intervals, to produce a plurality of 

20 carrier phase measurements and a plurality of carrier 

21 phase estimates at each of a succession of time points; 

22 determining the difference between the changes in 

23 carrier phase measurements and the changes in carrier 

24 phase estimates over successive time points for each 

25 satellite/receiver pair, to produce a plurality of 

26 error values for each time point; 



27 
28 
29 



defining each of the plurality of error values for each 
time point to be equal to predetermined function of the 
signal's sensitivity to small changes in the position 
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30 coordinates of the corresponding receiver and in the 

31 clocks of the corresponding receiver and satellite; 

32 forming and updating at each time point a predetermined 

33 matrix based on the plurality of error equations for 

34 each time point; 

35 the step of forming and updating including a step of 

36 eliminating dependence on previous values of the 

37 receiver and satellite clocks; and 

38 solving the matrix to determine the particular position 

39 coordinates for the first receiver that minimize the 

40 mean square value of the plurality of error values. 

1 2. A method as defined in claim 1, wherein: 

2 the step of forming includes a step of forming a square 

3 root information filter matrix using a triangular 

4 Householder algorithm; 

5 the matrix entries indicating sensitivity to the 

6 receiver and satellite clocks are located in a 

7 predetermined section of the matrix; and 

8 the step of updating includes a preliminary step of 

9 zeroing the predetermined section of the matrix. 

1 3. A method as defined in claim 1, wherein the step of 

2 forming includes a step of modifying the successive 

3 sets of error equations to make the sets substantially 

4 uncorrelated with each other. 



1 
2 
3 



4. 



A method as defined in claim 3 , wherein the step of 
modifying includes a step of whitening the successive 
sets of error equations. 
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1 5. A method as defined in claim 1 and further including 
steps following the step of solving, of: 



2 



3 
4 
5 
6 
7 

8 
•9 



11 
12 
13 

1 

2 
3 

1 
2 

3 
4 



6 
7 
8 
9 

10 
11 



modifying the plurality of error equations to make the 
successive equations for each satellite/receiver pair 
substantially uncorrelated with each other, thereby 
producing a plurality of modified error equations for 
each time point; 

forming and updating the predetermined matrix based on 
the plurality of modified error equations for each time 



10 point; and 



5 pair; and 



solving the matrix to determine the particular position 
coordinates for the first receiver that minimize the 
mean square value of the plurality of error values. 

A method as defined in claim 5, wherein the step of 
modifying includes a step of whitening the successive 
sets of error equations. 

A method as defined in claim 5, and further including 
steps, following the second step of solving, of: 

resolving any ambiguity in the integral number of 
carrier cycles between each satellite and receiver 



determining the particular position coordinates for the 
first receiver that minimize the mean square value of 
the plurality of error values, taking into account the 
integral number of carrier cycles determined, in the 
step of resolving, to be between each satellite and 
receiver pair. 
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L 8. A method as defined in claim 4, wherein the matrix used 

2 in the step of forming remains the same size following 

3 each updating. 

1 9. A method as defined in claim 1, wherein the matrix used 

2 in the step of forming remains the same size following 

3 each updating. 

I. 10. Apparatus for determining the position coordinates of a 

2 first predetermined site using signals broadcast from a 

3 plurality of orbiting satellites, the apparatus 

4 comprising: 

5 a first receiver located at the first site; 

6 a second receiver located at a second, reference site 

7 having known position coordinates; 

8 wherein the first and second receivers both have 

9 internal clocks and are adapted to receive the signals 

10 broadcast from the plurality of satellites and measure 

11 the carrier phase of each signal relative to their 

12 respective internal clocks, to produce a plurality of 

13 carrier phase measurements for each of a succession of 

14 time points; 

15 means for calculating the expected phase of each 

16 carrier signal received by the first and second 

17 receivers from each satellite, based on the satellite's 

18 orbit and an estimate of the receiver's position, to 

19 produce a plurality of carrier phase estimates for each 

20 of the succession of time points; 

21 matrix means for forming a predetermined matrix that 
.22 defines each of the plurality of error values to be 

23 equal to a predetermined function of its own 

24 sensitivity to small changes in the position 



WO 87/06410 



PCT/US87/00825 



25 

25 coordinates of the corresponding site and in the 

26 internal clocks of the corresponding receiver and 

27 satellite; 

28 means for determining changes in the carrier phase 

29 measurements and changes in the carrier phase estimates 

30 over successive time points, and for determining the 

31 differences between corresponding pairs of such' 

32 measurements and estimates for each satellite/receiver 

33 pair and for each time point, to produce a plurality of 

34 error values for each time point; 

35 wherein the matrix means includes means for updating 

36 the matrix to reflect the plurality of error equations 

37 for each time point, the updating means including means 
3 8 for eliminating dependence on previous values of the 

39 receiver and satellite clocks; and 

40 means for solving the matrix to determine the 

41 particular position coordinates for the first site that 

42 minimize the mean square value of the plurality of 

43 error values. 

1 11. Apparatus as defined in claim 10, wherein: 

2 the predetermined matrix formed by the matrix means is 

3 a square root information filter matrix; 

4 the matrix means includes means implementing a 

5 triangular Householder algorithm . to incorporate the 

6 successive error equations into the matrix; 

7 the matrix entries indicting sensitivity to the 

8 receiver and satellite clocks are located in a 

9 predetermined section of the matrix; and 
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1° the matrix means further includes means for zeroing the 

11 predetermined section of the matrix. 

1 12. Apparatus defined in claim 10, wherein the matrix means 

2 further includes means for modifying the successive 

3 sets of error equations to make the sets substantially 

4 uncorrelated with each other. 

1 13. Apparatus as defined in claim 12, wherein the means for 

2 modifying includes means for whitening the successive 

3 sets of error equations. 

1 14. Apparatus as defined in claim 10, and further 

2 including: 

3 means for modifying the plurality of error equations to 

4 make the successive equations for each 

5 satellite/receiver pair substantially uncorrelated with 

6 each other; thereby producing a plurality of modified 

7 error equations for each time point; 

8 means, operable after the first means for solving, for 

9 forming and updating a second predetermined matrix 

10 based on the plurality of modified error equations for 

11 each time point; and 

12 means for solving the second matrix to determine the 

13 particular position coordinates for the first receiver 

14 that minimize the mean square value of the plurality of 

15 error values. 

1 15. Apparatus as defined in claim 14, wherein the means for 

2 modifying includes means for whitening the successive 

3 sets of error equations. 



1 16. 

2 



Apparatus as defined in claim 14, and further 
including: 
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3 means, operable after the second means for solving, for 

4 resolving any ambiguity in the integral number of 

5 carrier cycles between each satellite and receiver 

6 pair; and 

7 means for determining the particular position 

8 coordinates for the first receiver that minimize the 

9 mean square value of the plurality of error values, 

10 taking into account the integral number of carrier 

11 cycles, determined by the means for resolving, to be 

12 between each satellite and receiver pair. 

1 17. Apparatus as defined in claim 14, wherein the first and 

2 second predetermined matrices remain the same size 

3 following each updating. 



1 
2 
3 



18. 



Apparatus as defined in claim 10, wherein the 
predetermined matrix remains the same size following 
each updating. 
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